Steady state velocity distributions of an oscillated granular gas 
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We use a three-dimensional molecular dynamics simulation to study the single particle distribution 
function of a dilute granular gas driven by a vertically oscillating plate at high accelerations (15g — 
90g). We find that the density and the temperature fields are essentially time-invariant above a 
height of about 35 particle diameters, where typically 20% of the grains are contained. These grains 
form the nonequilibrium steady state granular gas with a Knudsen number unity or greater. In the 
steady state region, the distribution function of horizontal velocities (scaled by the local horizontal 
temperature) is found to be nearly independent of height, even though the hydrodynamic fields 
vary with height. The high energy tails of the distribution functions are described by a stretched 
exponential ~ exp(— Be"), where a depends on the normal coefficient of restitution e (1.2 < a < 1.6), 
but q does not vary for a wide range of the friction parameter. We find that the distribution function 
of a frictionless inelastic hard sphere model can be made similar to that of a frictional model by 
adjusting e. However, there is no single value of e that mimics the frictional model over a range of 
heights. 

PACS numbers: 45.70.-n, 05.20.Dd, 05.70.Ln, 83.10.Rs 



I. INTRODUCTION 

A dilute gas in thermal equilibrium is sufficiently char- 
acterized by the pressure and temperature and is de- 
scribed by a simple relation, the equation of state. How- 
ever, when a gas is far from equilibrium, there is no gen- 
eral, finite set of variables specifying the state. The sin- 
gle particle distribution function /(r, v, t) is often suf- 
ficient to characterize the statistical properties of a di- 
lute nonequilibrium gas when correlations are negligible. 
Given this function, other quantities, such as moments of 
the distribution and transport coefficients, can be eval- 
uated. Dilute granular materials subject to an external 
forcing exhibit gaseous behaviors that share many analo- 
gies with a molecular gas, and they are often called gran- 
ular gases. Such a granular gas is always far from equilib- 
rium due to the dissipative collisions, and the deviation of 
/(r, v, t) from the Maxwell-Boltzmann (MB) distribution 
has been of great interest in recent years 0, 0, 0, 0, Q ■ 

The velocity distributions of a vibro-fluidized granular 
gas were first measured by Warr et al. they stud- 
ied the distribution functions of grains confined between 
two transparent plates and concluded that the distribu- 
tion was consistent with the MB distribution function. 
Recently, the same system has been studied by Rouyer 
et al. 5J, who found a universal distribution function 
of the form ~ exp(— £>|v| 15 ) for the entire range of ve- 
locities studied, where B was a parameter. The authors 
reported that this functional form fit their measurements 
for a wide range of oscillation parameters for various ma- 
terials; thus the granular temperature, the second mo- 
ment of the distribution, was the only parameter of the 
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distribution function. 

There have been numerical studies of vibrated inelas- 
tic hard disks, subject to a saw-tooth type oscillation, in 
the presence of gravity @ and in the absence of grav- 
ity 0| . Such forcing is often used in theoretical studies 
as a simplification of the sinusoidal oscillation used in 
experiments, assuming that the asymptotic behavior of 
the hydrodynamic fields far from the plate is the same; 
however, it is not known a priori how far from the oscil- 
lating plate one must be in order for this assumption to 
be valid. 

In this paper, we perform a simulation that is as close 
as possible to three-dimensional experiments on verti- 
cally oscillated granular gases. We use a previously val- 
idated molecular dynamics (MD) simulation 8, 9]. The 
hydrodynamic fields are oscillatory near the plate, and 
their oscillatory behavior decays with height. Above 
some height, the fields are not correlated with the os- 
cillation of the plate and are essentially time-invariant. 
We study the distribution functions in this nonequilib- 
rium steady state region. To focus on the distributions 
due to the intrinsic dynamics of the granular gas, we do 
not impose sidewalls or include air. We also study how 
the distribution changes with the friction. In many the- 
oretical or numerical studies of granular fluids, granular 
materials are modeled as frictionless inelastic hard disks 
or spheres; however, no granular materials are friction- 
less, in the same way that none of them are elastic. We 
check if the role of friction can be incorporated into the 
inelasticity by adjusting the value of the normal coef- 
ficient of restitution. In this paper we discuss the dis- 
tributions only in the steady state region; those in the 
oscillatory region near the plate will be discussed in a 
separate paper [lj} . 

The rest of the paper is organized as follows. In Sec- 
tion II, the system under consideration, the data analysis 
method, and the collision model are described. Results 
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are presented in Section III and discussed in Section IV. 



B. collision model 



II. METHOD 
A. System and data analysis 

We use both a frictional and frictionless, inelastic hard 
sphere MD simulation. We consider 133 328 monodis- 
perse spheres of unit mass and of diameter a = 165 \xm 
in a container with square bottom of area 200cr x 200cr 
(the average depth of the layer at rest is approximately 
3er), where periodic boundary conditions are imposed in 
both horizontal directions. We choose the same particle 
size as in Ref . , as the patterns were quantitatively re- 
produced for a wide range of oscillation parameters with 
this particle size; however, as long as the collision model 
is valid, all the length scales can be normalized by a. We 
assume the bottom plate of the container is made of the 
same material as grains; we use the same material coeffi- 
cients for the inelasticity and the friction as grains. The 
bottom plate is subject to a vertical sinusoidal oscilla- 
tion with an amplitude A and a frequency /. We vary 
the oscillation parameters in the range of 3<7 < A < 10a 
and 40 Hz < / < 170 Hz, which approximately corre- 
sponds to 0.35 m/s < V m ax (= 2irfA) < 0.75 m/s and 
15.g < a max [= A(2nf) 2 ] < 90g, where g is the accelera- 
tion due to gravity. We check that with our parameters 
no mean flow develops and that grains rarely reach to the 
top, which is fixed at 300<r. 

Hydrodynamic fields and the distribution functions are 
analyzed by binning the box into horizontal slabs of 
height a, as the system is invariant under the transla- 
tion in both horizontal directions, in the absence of any 
mean flow. We use the granular volume fraction v for 
the density, which is the ratio of the volume occupied by 
grains to the volume of each horizontal slab. We consider 
the following three granular temperatures separately: 
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where x and y are horizontal directions that are indistin- 
guishable, z is the vertical direction, v is a velocity vector 
for each grain, and the ensemble average ( ) is taken over 
the particles in the same bin at the same phase angle dur- 
ing 40 cycles, after initial transients have decayed. We 
define the scaled horizontal velocity to be 



We implement the collision model that was originally 
proposed by Maw et al. simplified by Walton |12|, 
and then experimentally tested by Foerster et al. |13|. 
This model updates the velocity after a collision accord- 
ing to the three parameters, the normal coefficient of 
restitution e (e [0, 1]), the coefficient of friction [i, which 
relates the tangential force to the normal force at colli- 
sion using Coulomb's law and then determines the tan- 
gential coefficient of restitution (3 (e [—1,1]), and the 
maximum tangential coefficient of restitution /3o, which 
represents the tangential restitution of the surface veloc- 
ity when the colliding particles slide discontinuously at 
the contact point. 

At collision, it is convenient to decompose the relative 
colliding velocities into the components normal (v ra ) and 
tangential (v t ) to the normalized relative displacement 
vector fi2 = (ri — r 2 )/|i"i — r 2 |, where ri and r 2 are 
displacement vectors of grains 1 and 2, and the same 
notation is used for v: 



V„ = (V12 • fi 2 )fi2 = V n ii2, 
V t = ?12 X (V12 X ?1 2 ) = V12 



(6) 
(7) 



The relative surface velocity at collision, v s , for monodis- 
perse spheres of diameter a is 



-fi2 X (wi + W 2 ) = V S V S 



(8) 



where wi and W2 are the angular velocities of grains 1 
and 2, respectively. 

For monodisperse spheres of diameter a and unit 
mass, the linear and angular momenta conservations and 
the definitions of the normal coefficient of restitution 
e = —v*Jv n and the tangential coefficient of restitution 
(3 = —v*/v s (post-collisional velocities are indicated by 
superscript *, and pre-collisional values have no super- 
script) give the changes in the velocities at the collision: 



Avi„ = 
Av u = 
Awi = 



1 / 

-Av 2 „ = - (1 + e) v„, 
K(l + 0) 

' Av2t = 2(kTT) Vs ' 

(1+/?) 
<t(K + 1) 



-A\V2 



ri2 x v s 



(9) 
(10) 

(11) 



where K = 41 /a 2 is a geometrical factor relating the mo- 
mentum transfer from the translational degrees of free- 
dom to rotational degrees of freedom, and / is the mo- 
ment of inertia about the center of the grain. For a uni- 
form density sphere, K is 2/5. 

We use a velocity-dependent normal coefficient of resti- 
tution as in Ref. 8], to account for the viscoelasticity of 
the real grains: 



= K - {v x ))h 
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where eo is a positive constant less than unity. Since 
we impose high forcing (V max > 0.35 m/s while ^fga = 
0.04 m/s) the collision probability for relative colliding 
velocities v n less than ^fga is small, and using a velocity- 
independent e does not result in any noticeable difference, 
compared to using a velocity-independent one e = eo ; the 
same was true in Ref. [l^J. We use the symbol e for eo 
hereafter. 

In collisions of real granular materials, not only is the 
relative surface velocity reduced, but also the stored tan- 
gential strain energy in the contact region can often re- 
verse the direction of the relative surface velocity. To 
account for this effect, the tangential coefficient of resti- 
tution (3 could be positive, leading to the range of 
as [—1,1]. There are two kinds of frictional interaction 
at collisions, sliding and rolling friction, which are ac- 
counted for by the following formula for 0: 







(13) 



where 0o is the maximum tangential coefficient of resti- 
tution. For sliding friction, the tangential impulse is as- 
sumed to be the normal impulse multiplied by /i. When 
is identically negative unity (or simply \i = 0), this model 
reduces to the frictionless interaction. For the special 
case v s = 0, the collision is treated as frictionless. This 
friction model is still a simplification of the real frictional 
interaction; there is no clear-cut distinction between the 
two types of frictions for real grains, and even a transfer 
of energy from the rotational to translational degrees of 
freedom, which results in e larger than unity, has been 
observed However, this collision model is accurate 
enough to reproduce many phenomena, including stand- 
ing wave pattern formation in vertically oscillated gran- 
ular layers, when the parameters are properly chosen. In 
Refs. |j| and Q, e = 0.7, O = 0.35, and fi = 0.5 were 
used. 



III. RESULTS 

A. Hydrodynamic fields and steady state 

Due to the oscillatory boundary forcing, the hydrody- 
namic fields, the volume fraction v and the granular tem- 
peratures (T, T x , and T z ), depend on height z and time 
t (Fig. Q near the oscillating plate; the temperatures 
exhibit stronger oscillatory behaviors than the density. 
Since the energy is injected mainly through the verti- 
cal velocities, the granular temperature is anisotropic, as 
illustrated in Fig.Q T z is larger than its horizontal coun- 
terpart T x , and the former is significantly larger near the 
bottom plate, where the hydrodynamic fields are oscilla- 
tory. The vertical temperature increases almost linearly 
with height for zj a > 35; however, the temperature T in- 
creases slower than linearly, as the slope of T x decreases 
with height and T x levels off for zja > 120 (Fig. QJ. 
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FIG. 1: The volume fraction v (thick solid line) and the 
granular temperature T (thick dashed line) as a function of 
height at different times during a cycle, where e = 0.9, (3o = 
0.35, fx = 0.5, Vmax = 0.55 m/s (a ma x = 60#, A = 3cr, 
and / = 169 Hz), and ft is set to zero when the plate is at 
the equilibrium position moving upward. Above some height 
(2/17 ~ 35), the hydrodynamic fields do not vary much in 
time. The horizontal temperature T x (thin dashed line) is 
smaller than the vertical temperature T z (thin solid line). The 
container bottom is indicated by the vertical gray line. 



A similar increase of the temperature with height was 
observed in an open system of frictionless inelastic hard 
disks or spheres subject to a thermal bottom heating 0] 
and a saw-tooth type vibration 0] . We characterize the 
oscillation parameters only by V max , as we observe for the 
parameters in our study that the hydrodynamic fields in 
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FIG. 2: This superposition of the volume fraction and tem- 
perature fields at five different times in a cycle (see Fig. 
illustrates that the fields are nearly time-independent above 
z/a ~ 35, which we call the steady state region. 



FIG. 4: The mean free path A, estimated from a formula for a 
gas of hard spheres, and the Knudsen number Kn, estimated 
by using A and the length scale of the density l v (see text), in 
the steady state region. 



the steady state are nearly the same for the same V max , 
even for different combinations of a m ax and /; such scal- 
ing behavior was also observed in Ref . [l^ ■ 

During each cycle, a normal shock forms at the im- 
pact from the bottom plate and propagates upward 0] . 
As the shock propagates, it decays and becomes unde- 
tectable above some height (z/a « 35), rather than prop- 
agating up through the entire granular media (which was 
the case in Ref. |l8|). Above this height, the hydrody- 
namic fields are invariant in time, and the granular gas 
forms a nonequilibrium steady state (Fig.[2J|, where about 
20% of the grains are contained in this case; this fraction 
depends on the oscillation parameters. 

With the parameters used in this paper the granular 
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FIG. 3: The number of grain-grain collisions per grain during 
a cycle N co u (solid line), and time-averaged volume fraction 
Vavg (dot-dashed line, multiplied by 100) and time-averaged 
granular temperature T avg (dashed line, multiplied by 100), 
over sixty different equally spaced times during a cycle. N co u 
is less than 1 in the steady state region (z/a > 35). Inset : 
The same quantities (without multiplications) on a logarith- 
mic scale. 



temperatures are nonzero throughout the cycle, as grains 
do not solidify after the shock passes through, in contrast 
to the case in Ref. [3. As a result, when the bottom 
plate moves down, the granular gas expands, and an ex- 
pansion wave propagates downward (see the temperature 
peaks near the plate for ft > 0.4 in Fig. 

We count the number of grain-grain collisions per grain 
during a cycle ( N co u in Fig. 0) , and find that N co u is less 
than unity in the steady state region; the granular gas 
in the steady state is nearly collisionless. We estimate 
the mean free path using the formula for a gas of hard 
spheres \{z)/a = ^V^naf)" 1 = y/n '/ (l2-y/2v) (where 
n is the number density) which ranges between 5 
and 280 for 40 < z/a < 100 (Fig. gj). When we esti- 
mate the mean free path using the measured collision 
frequency and the thermal speed, we get a similar result. 
We fit the density with piecewise exponential functions 
[~ exp(— [z — Zo)/l v )] in the steady state region, and ob- 
tain a hydrodynamic length scale l u /a between 12 and 
15 for 40 < z/a < 100. We obtain a similar length scale 
from piecewise linear fitting of the temperature T in the 
same region. We calculate the Knudsen number Kn, de- 
fined as the ratio of the mean free path to the length scale 
of the macroscopic gradients [20j. using X(z) and l v \ Kn 
ranges from 0.5 to 20 in the region 40 < z/a < 100 
(Fig. II). 



B. Height-independence of the distribution 

The distribution of scaled horizontal velocities should 
be symmetric as a consequence of the symmetry of the 
system. We calculate the skewness of the distribution, 
73 = M.'i/M^' 1 ' , where M n is the n th moment of the 
distribution 



M n = / c™f(c x )dc x , 
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FIG. 5: Above some height (z/a ~ 42, indicated by a ver- 
tical dashed line), the kurtosis 74 is nearly time-invariant. 
The horizontal granular temperature T x (dashed lines) and 
74 (solid lines) of the horizontal velocity distribution function 
are shown at five different times during a cycle (cf. Fig.0. 



FIG. 7: Kurtosis for three values of e, as a function of height 
at ft = —0.2, where fj, and /3a are set to 0.5 and 0.35. Different 
forcings are applied for each case to achieve similar profiles 
of the hydrodynamic fields: V m ax (a m ax) = 0.4 m/s (43<?), 
0.55 m/s (6O5), and 0.66 m/s (72g) for e = 0.95, e = 0.9, and 
e = 0.85, respectively. 



and check that |7 3 | is less than 0.01 for all the dis- 
tributions we study. The lowest order deviation from 
the MB distribution is characterized by the flatness of 
the distribution, which is called the fourth cumulant 
or the kurtosis. It was used to quantify the devia- 
tion from the MB distribution of the homogeneously 
cooling state [2TL l22l l23l |. the homogeneously heated 
state [2l|, I24T and granular gases subject to a bound- 
ary forcing |oll25|. The kurtosis 74 (= M4/M?, - 3) is 
defined so that it vanishes for the MB distribution, and 
we find that it also does not change in time above some 
height (Fig. [SJ. Further, in the steady state region, 74 is 
nearly independent of the height, even though both the 
density and the temperature change; the distributions 
at different heights in the steady state region are hardly 
distinguishable (Fig. [SJ. Also, the kurtosis in the steady 
state region does not vary much for a wide range of the 
oscillation parameters: for 0.35 m/s < V rnax < 0.75 m/s 
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FIG. 6: The horizontal velocity distribution functions at four 
different heights (compared with /ms , the solid line) obtained 
at ft = —0.2. There is no noticeable difference among the 
distributions in the steady state region, 40 < z/a < 80. 



(and other parameters fixed), 74 changes less than 10%. 
A similar absence of height dependence of the velocity 
distribution function was found in a recent ex per iment 
on a vertically oscillated quasi-2d granular gas [2q . 

C. Velocity distributions 

We first examine the dependence of the distribution 
on e; we measure 74 for three different values of e, while 
[3q and \i are kept at 0.35 and 0.5, respectively. We find 
that 74 significantly decreases with increasing e in the 
oscillatory state, however, it decreases only slightly in 
the steady state region (Fig. 01 . 

Now we compare our results with the MB distribution 
of variance l/\/2, 

!mb{c x ) = -=exp(-cl), (15) 

The steady state distributions obtained for the parame- 
ters in Fig. are overpopulated in the high energy tails 
and underpopulated at small velocities (Fig. |SJ|, com- 
pared to fMB- The differences of the distributions for 
various e's are hardly noticeable both on linear and loga- 
rithmic scales (Fig. [8}, but on a double logarithmic scale 
plot the tails of the distributions are described by dif- 
ferent functions (Fig. EJ. We investigate the functional 
form of the distributions by fitting them (after the nor- 
malization by the value at c x = 0) with a stretched ex- 
ponential function exp(— c"). We find that the exponent 
a changes from 1.9 (indicated by dashed lines) to some 
smaller value (solid lines), depending on e, as the velocity 
increases. We have not investigated lower values of e to 
avoid issues such as cluster formation. 

We now keep e and (3q at 0.9 and 0.35, respectively, 
and change the value of [i. The profile of 74 in the os- 
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FIG. 8: The distribution functions of scaled horizontal veloc- 
ities f(c x ) for the cases in Fig.0 on linear (top panel) and log- 
arithmic (bottom panel) scales. Although the kurtosis slightly 
decreases with e, the difference between the distributions is 
hardly distinguishable on both scales. The ranges for the av- 
eraging in height were between 30ct and 45a for e = 0.95, 40a 
and 55a for e = 0.9, and 45a and 60a for e = 0.85; the averag- 
ing was done over relatively similar heights in the steady state 
regions (see Fig. 0. The solid line is the MB distribution. 
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FIG. 9: Double logarithm of the distributions of scaled hor- 
izontal velocities for the cases in Fig. as a function of the 
logarithm of c x . The slope corresponds to the negative expo- 
nent, —a, of a stretched exponential function exp(— c"). a is 
the same for small velocities for three different e's, however, 
it depends on e in high energy tails. Dashed lines correspond 
to a — 1.9, and the solid lines (from the top) correspond to a 
= 1.2, 1.4, and 1.6, respectively (indicated by the numbers). 



from the frictional model. 

The rotational kinetic energy is two orders of magni- 
tude smaller than its translational counterpart for the 
cases studied in this paper. However, the presence of the 
friction reduces the expansion of the granular gas signif- 
icantly, because the friction reduces the mobility of the 
grains and increases the collision frequency [271 ]. The 
mean height of frictional inelastic hard spheres exhibits 
a different scaling behavior with the plate velocity from 
that of frictionless spheres. Only the frictional sphere 
model reproduces the experimental observations p7l |28| . 

The 74 's obtained from the simulations of frictionless 



cillatory region changes significantly with /i, however, it 
is nearly unchanged in the steady state region (Fig. llOjl . 
The velocity distributions in this region for three differ- 
ent values of \i in Fig. ^] are also hardly distinguishable. 
We observe that the distribution function depends also 
on the density, as in Refs. 0,013; however, we do not 
investigate this dependency systematically. 



D. Frictionless inelastic hard sphere model 

In theoretical and numerical studies, granular mate- 
rials are often modeled as smooth (frictionless) inelas- 
tic hard disks or spheres, assuming that the friction is 
a secondary effect that can be neglected or that both 
the inelasticity and the friction can be incorporated to- 
gether into a so-called effective coefficient of restitution. 
In this Section, we discuss how the velocity distribution 
changes when the friction is not included, and we show 
that the frictionless model exhibits qualitative differences 
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FIG. 10: Kurtosis of the distributions of scaled horizontal 
velocities as a function of height at ft = —0.2, for three values 
of fi (/3o and e are kept at 0.35 and 0.9, respectively). In the 
oscillatory region, 74 increases with fi; however, 74 is nearly 
the same within the uncertainty in the steady state region. 
The same forcing (V ma x = 0.55 m/s) is applied to the three 
cases. 
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FIG. 11: Kurtosis of the distributions of scaled horizontal 
velocities of frictionless spheres, as a function of height for 
three values of e obtained at ft = —0.2. Different forcings (the 
same as in Fig. [7J are applied for each case: V max (a-max) = 
0.4 m/s (43s), 0.55 m/s (60g), and 0.66 m/s (72g) for e = 0.95, 
e = 0.9, and e = 0.85, respectively. 



particles for the same forcings as in Fig. [7| are illustrated 
in Fig. ^] In the absence of friction, the layer expands 
much more, and the steady state occurs at greater height, 
z/a > 100. In both the oscillatory and the steady state 
regions, values of 74 are smaller than those of frictional 
spheres (compare Fig. ^2 with Figs. [3 and llOfl ; the dis- 
tribution deviates from /mb only slightly. The kurtosis 
decreases with increasing e, and the difference among the 
distributions for the three e's are small. These distribu- 
tions have four crossovers with Jmb', they are overpopu- 
lated both at very small and high velocities and are un- 
derpopulated in between, compared to /m b ■ We fit them 
with a stretched exponential function, and find that a is 
2.0 for small velocities, and that it depends on e for the 
high energy tails (Fig. I12fl . as in frictional hard spheres. 

Since the functional form of the distribution depends 
on e, we can get a similar distribution function for the 
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FIG. 12: Double logarithm of the rescaled horizontal velocity 
distribution functions (normalized by its value at zero) for the 
cases in Fig. 1111 as a function of the logarithm of c x . Dashed 
lines correspond to a — 2.0, and the solid lines (from the top) 
correspond to a = 1.65, 1.85, and 1.9, respectively. 
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FIG. 13: The volume fraction v (gray scale and contour lines) 
at ft = 0.25 as a function of height and v z [(a) and (c)], and 
of height and v x [(b) and (d)], obtained from simulations of 
frictional hard spheres [(a) and (b); e = 0.9, n = 0.5], and 
of frictionless hard spheres [(c) and (d); e = 0.7, fi = 0] 
at Vmax = 0.55 m/s; e is adjusted in the frictionless case 
to obtain comparable overall dissipation and similar velocity 
distribution functions in the steady state region. The friction- 
less spheres spread more smoothly in height, and they do not 
yield the sharp gradient in the density around z/a ~ 13 as in 
the frictional case, no matter what the value of e is; the den- 
sity profile of the frictionless spheres is qualitatively different 
from that of frictional spheres. Contour lines correspond to 
(0.03,0.3,0.6,0.9,1.2,1.5) xlO" 3 from outside. 



steady state by adjusting e. For instance, for jj, = 0, 
e = 0.7, and V max — 0.55 m/s (the same forcing as in 
Fig. EJ, we obtain 74 w 0.5 for the steady state region; 
the steady state distribution is similar to the one in Fig.0 
for e = 0.9 and /1 = 0.5. However, we find that no single 
value of e mimics the hydrodynamic fields or the distri- 
bution function of frictional hard spheres, both in the 
oscillatory and steady state regions; the effect of friction 
cannot be taken over by an adjusted normal coefficient of 
restitution. The difference between the results of the two 
models is illustrated in Fig. ^3 where velocities are not 
rescaled for better comparison. The outermost contour 
lines in both models become similar when e in the fric- 
tionless model is adjusted as a free parameter [compare 
Figs. I13t a) and (c), or (b) and (d)]; however, the overall 
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shape of the density contours cannot be matched by ad- 
justing only e. Note that the density changes rapidly with 
height and c>1.5x 10 -3 at z/a w 10 near v z « v x « 
in the frictional model, whereas in the frictionless model, 
the particles spread more smoothly in height, and there 
is no region for v > 1.2 x 1CP 3 . 

IV. CONCLUSIONS 

We have studied the horizontal velocity distribution 
function of vertically oscillated dilute granular gas, using 
a molecular dynamics simulation of frictional, inelastic 
hard spheres. The hydrodynamic fields are oscillatory 
in time near the oscillating bottom plate due to a shock 
wave and an expansion wave. However, the helds are 
nearly stationary above some height, thus constituting 
a granular gas in a nonequilibrium steady state. The 
steady state region forms a granular analog of a nearly 
collisionless Knudsen gas (Figs. and @J. We find that 
the dependence of the distribution functions in this gran- 
ular Knudsen gas regime on the forcing and material pa- 
rameters is very weak, even though the distributions in 
the collisional bulk at lower heights depend strongly on 
the forcing and material parameters (Figs. 171 and HUt . 
The behavior of an ordi nary Knudsen gas is determined 
by boundary conditions |20|. Although we do not know 
whether boundary conditions or collisions are dominant 
in determining the behavior of our granular Knudsen gas, 
we note that this gas does not depend much on the prop- 
erties of its only boundary, which is the oscillatory region 
close to the plate. 

The functional form of the horizontal velocity distri- 
bution in the steady state region is nearly independent 
of height, when velocities are scaled by horizontal tem- 
perature (Fig. |HJ), even though the hydrodynamic fields 
continue to change. The distribution function is broader 
than the MB distribution, being underpopulated at small 
velocities and overpopulated in the high energy tails 
(Fig. |SJ . We do not observe a universal functional form 
for the distribution function (Fig. |5J). The functional 
form of the high energy tail changes with the dissipa- 



tion parameters (e and fi) and the oscillation parameter 
(V max ). The dependence on \x in the steady state region 
is very weak (Fig. 110(1 . 

Our conclusions regarding the absence of a universal 
distribution function differ from that of Ref. [5|, because: 
(1) We studied the local distribution function, while in 
Ref. [f| the authors obtained the distribution by averag- 
ing over space and time; they assumed that the spatial 
and temporal variation was negligible near the center of 
the oscillating box, based on their observation of a weak 
dependence of the density. We find that the time depen- 
dence of the density is weak, but that of the temperature 
is strong in the oscillatory state region (Figs. and [21 . 
Note that if a distribution is averaged over different tem- 
peratures, even the MB distribution function leads to a 
different resultant distribution function. (2) Our system 
is different from that in Ref. we do not have either 
air or sidewalls, and our container is much taller, so that 
the bottom plate is the only energy source in our case. 
How air and sidewalls affect the dynamics of a granular 
gas is yet to be clarified. 

We also studied the velocity distributions of friction- 
less inelastic hard spheres, and examined the possibility 
of including frictional effects using an effective normal 
coefficient of restitution. We found that no single effec- 
tive restitution coefficient could describe the frictionless 
gas at different heights. 

Velocities of a granular gas, even in the dilute limit, are 
strongly correlated, and the correlations depend on the 
density and the coefficient of restitution [24| . The depen- 
dence of the distribution on the density implies that the 
single particle distribution of a dilute granular gas can- 
not play a role equivalent to that in a dilute ordinary gas; 
it is not sufficient to specify statistical properties of the 
gas. However, the knowledge of the single particle distri- 
bution of this complex nonequilibrium gas is still of great 
importance for the purpose of the first approximation. 
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